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ABSTRACT 

In the ellipsoidal collapse model, the critical density for the collapse of a gravitationally 
bound object is a function of its mass. In the excursion set formalism, this translates 
into a moving barrier problem such that the mass function of dark matter haloes 
is given by the first crossing probability of a random walker across the barrier. In 
this paper, we study this first crossing probability analytically. Complete solutions 
are obtained for barriers that vary as square root of time and square of time. Large 
and small time asymptotic behaviour is derived. For arbitrary power-law barriers, the 
large time behaviour is determined. The solutions allow us to derive useful inferences 
about the scaling of the conditional mass function in terms of present day halo masses 
and look back redshifts. As an application of our results, we compare the estimates 
of major merger rates of haloes in constant and moving barrier models and find that 
for massive haloes (10 12_13 Mq) the latter predicts significantly higher merger rates 
towards high redshifts [z > 4). 

Key words: galaxies: clusters: general - cosmology: theory - dark matter. 



1 INTRODUCTION 

The problem of determining the mass function of gravita- 
tionally bo und structures was first a ddressed in a successful 
manner bv lPress fc Schechterl <ll974f) . whose model assumed 
that the primordial density fluctuations filtered on a given 
mass scale were Gaussian distributed. Since their model, 
many detailed schemes have been investigat ed and per- 
fected . The excursion set appro ach developed in Bond et. alJ 
Jl99lD (hereafter BCEK) and lLacev fc Colel (Il993ft (here- 
after LC93) is a convenient tool for deriving the uncondi- 
tional and the conditional mass functions within the frame- 
work of Gaussian random fields. The approach has been suc- 
cessfully used to cre ate Monte-Carlo merging history trees of 
dark matter haloes (Kauffmann. White & Guide rdonilll993r 
ISomerville fc Kolattlll999USheth fc Lemsonlll999T) . As com- 
pared to TV -body simulations, Monte-Carlo merger trees 
provide alternative fast er methods fo r study ing the build-up 
of dark matter haloes. iMo fc White! 1(1996) further showed 
how the spatial distribution of haloes may be quantified 
within the excursion set approach. 

The excursion set approach is based on the following 
principles. Consider a dark matter inhomogeneity centred 
around some point in the universe. The smoothed density 
contrast within a radius R around this point is defined as 
S(R) = [p(R) — po]/po, where p(R) is the density of mat- 
ter within R and po is the mean background density of the 
universe. If the density contrast is greater than the crit- 
ical density for collapse, the matter contained within the 



volume eventually collapses to form a bound object. Practi- 
cally, S(R) is obtained by convolving the matter density field 
with some spherically symmetric function Wn(r) of radial 
extent R. The variance of the smoothed density contrast is 
then (e.g., LC93) 



S(R) = (\S(R)\ 2 ) = 



i r 

Wo 



dk (\S k \ 2 )WUk) fc 2 , (i) 



where Sk's are the Fourier amplitudes of the field and Wn(k) 
is the Fourier transform of the window function Wr(t). 

A convenient choice for the window function is the 
sharp fc space function, where the cutoff wavenumber k 3 
is related to the Lagrangian radius (R) of an object by 
R ~ (97r/2) 1/3 fc s _1 (LC93). The equivalent mass scale is 
M ~ poR 3 - The sharp fc window function has the advantage 
that 8[R(k s )} executes a random walk with every increment 
in the size of the window. The linearly extrapolated critical 
density for collapse now serves as an absorbing barrier for 
random walk trajectories and the mass function of collapsed 
objects is given by the first crossing distribution of these 
random walks across the barrier. 

In the spherical collapse model, the critical density 
is independent of the collapsing halo masses. When lin- 
early extrapolated to the present day, it has a value of 
5 C = 1.686 in a standard Cold Dark Matter (CDM) uni- 
verse. The mass function of dark matter haloes derived us- 
ing the excursion set approach with spherical collapse model 
lie wi thin 10 — 30 % of t he results from N -body simula- 
tions dJenkins et al]|200lT) . Detailed comparison of simula- 
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tions with analytical results s hows discrepancies for both 
smal l and large mass haloe s <Gelb fe Bertschingerl Il994t 
lLacev fc Coldll994l; lTOTmen|ll998l : ISheth fe Tormedll999l) . 

ISheth. Mo fc Tormenl (1200 J) (hereafter SMT01) and 
iMonacoT lll997dlbD investigated a non-spherical alternative 
for collapse of over-densities. In particular they focused on 
an ellipsoidal collapse scenario which derives support from 
the tr iaxial nature of pe r turbations in Gaussi an density 
fields l lDoroshkevicrJll970l: Ifjardeen et. al.lll986fl . They ar- 
gued that the main effect of an ellipsoidal collapse is to 
introduce a dependence of the critical density on the halo 
mass. In the excursion set formalism this amounts to incor- 
porating a 'moving' barrier instead of a fixed one. 

In a subsequent paper ISheth fc Tormenl 1200 2) (here- 
after ST02) presented a detailed discussion on moving bar- 
rier models. They considered random walks in one dimen- 
sions x [= S(k s )] diffusing with time t [= S(k 3 )] starting at 
x — when t = 0. The first crossings with a moving bar- 
rier of the form B(t) = a + bt 1 were studied. Based on first 
crossing distributions f(t) obtained from Monte-Carlo sim- 
ulations, they suggested that for a barrier B(t), f(t) has the 
form 



/(*) = 



-S(t) 2 /(2t) 

2rf/ 2 



n. ^ 



(-1)" 



>n!( 7 - 



IJ(7- ^ 



(2) 

where n* ~ 5. In this expression, the diffusion constant D 
has been set equal to 1. Equation [2] when specialised to 
constant (6 = 0) and linear barriers (7 = 1), gives the correct 
answer. For these two barriers, the first passage distribution 
is easily obtained by the reflection principle felled Il97fl) 
and is 



fit) 




2Dt 

-(a + bt) 



0, 



S£] 7=1- 



(3) 



However, the validity of equation [5] for other kinds of barri- 
ers remains unchecked. For the two barrier problem, corre- 
sponding to the conditional mass function (the conditional 
mass function gives the progenitor mass distribution for a 
given present day halo at a given look-back redshift), fur- 
ther formulae were suggested based on a generalisation of 
the above expression. 

The precise value of 7 applicable to the ellipsoidal col- 
lapse model seems to lie between 0.5 and 1.0. Based on nu- 
merics and other arguments, SMT01 argue that 7 w 0.6 
with 6 > 0. On the other hand, by applying Lagrangian 
pertur bation theory and considering an ellipsoidal collapse 
model. iMonacol Jl997albl) concludes that 7 = 1/2 with b < 0. 
Also, from the study of tidal torques on galaxy evolution 
iDel PopoIo fc Gamberd dl998|) showed tha t barrier is of the 
form 6 > and 7 » 0.58 JPel Popoldl2002l) . 

In this paper, we present an analysis of the first passage 
distribution f(t) for a random walker with a moving barrier. 
In the process we test the validity of equation |2] We present 
an analytical solution of the square root barrier (7 = 1/2). 
This barrier is close to the one studied in SMT01 (7 w 0.6). 
We also solve for the quadratic boundary (7 = 2). The large 
and small time behaviour of these solutions are also derived. 
For arbitrary 7, we calculate the large time behaviour of the 
first crossing probability. Our results show that equation^ 
is not correct in general and needs to be modified. 



We also present the methodology for approaching the 
problem of conditional mass functions in the context of mov- 
ing barrier models. The conditional mass function requires 
solving a two barrier problem. Section |5] describes how to 
modify the moving barrier problem so that both the con- 
ditional as well as unconditional mass distributions may be 
obtained. 

The rest of the paper is organised as follows. In sec- 
tional we present the analytical solution for the square root 
barrier (7 = 1/2). The large and small time asymptotic 
behaviour is derived. Expressions are also obtained for con- 
ditional mass distributions. The analytical results are com- 
pared with results from Monte-Carlo simulations. 

In section 01 we consider another solvable limit - the 
quadratic barrier (7 = 2). The full solution and the large 
time asymptotic behaviour are derived, and compared with 
results from simulations. For both these cases, we show that 
the results are not consistent with the formula in equation|5| 

In section we consider barriers with arbitrary 7. For 
this general case, we argue what the large time behaviour 
of f(t) should be. Using special algorithms, we numerically 
obtain f(t) for large times and confirm our prediction. 

In section HJ we fit the halo mass distribution obtained 
from iV-body simulations with the square root barrier re- 
sults. While a good fit is obtained, the numerical data is not 
good enough to differentiate between different values of 7. 

In section|7] we describe an application of the results to 
estimate major merger rates of haloes. These rates are cal- 
culated using the conditional mass functions for the square 
root barrier. These rates are compared with those derived 
from a constant barrier Press- Schechter model. The analy- 
sis suggests that while the two models are very similar in 
describing the low redshift evolution of these rates, there 
is a systematic deviation to wards high redshifts. In partic- 
ular the cumulative major merger rates of massive haloes 
(10 12_13 M©) differ significantly at redshifts z > 4. 

We conclude with a summary and discussion in Sec- 
tion |H| In Appendix^ results for parabolic cylinder func- 
tions that are relevant to the paper are reproduced. 

Unless otherwise stated, the cosmological parameters 
used will be Q m = 0.3, £Ia = 0.7, h = 0.7 and root mean 
square fluctuations at 8 ft _1 Mpc as as = 0.9. We also inter- 
changeably use x <-> S and t «-» S. 



2 THE MOVING BARRIER MODEL 

For the ellipsoidal collapse model, the critical density con- 
trast depends on the varianc e or 'scale' S and has the form 
JSheth. Mo fc TormenlEoOlh 



5 ec (S, z) = y/q 5 sc (z) 



(4) 



Here, S„ c (z) is the spherical collapse critical density at red- 
shift z and St, = 6g C (z). The best-fit values of parameters /?, 
7 and q for the mass fu nction of haloes in G IF simulations of 
the Virgo consortium lljenkins et al.lfl998l) are 0.5, 0.6 and 
0.707 respectively. The value of q depends on the way haloes 
are identified in the simulations and therefore is a function 
of the link-length in the halo finding algorithms. 

The variable 8 executes a random walk when S is in- 
creased. We are interested in the probability that 5 exceeds 
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Figure 1. Illustration of the two barrier problem corresponding 
to the conditional mass distribution. 



the critical density S ec at scale S. In the following analysis 
it is convenient to substitute the variables S and S by vari- 
ables x and t respectively. Consider a random walk starting 
at x = when t = 0. Let P(x, t) be the probability of finding 
it at x at time t. Then, P(x,t) obeys the diffusion equation 



dP(x,t) _ Dd 2 P(x,t) 



(•») 



dt 2 dx 2 ' 

where D is the diffusion constant with a value one in the 
present context. In terms of x and t, the barrier in equation 
2]takes the form a + bt 1 , where a, b are constants. 

In what follows, we will consider modified barriers of 
the form 



B(t) =a + b(t + toT 



(6) 



where to is some constant (see below for motivation for this 
choice). Let g(a,b,to;t) be the probability that the random 
walk crosses the barrier B(t) for the first time at time t. 
The unconditional mass distribution is obtained by setting 
to = 0. To calculate the conditional mass function one has 
to consider a two barrier problem. 

Consider the illustration in figure The random walk 
trajectories begin at position O'. The moving barrier B\ 
represents the critical density at present epoch zi and has 
the form B\ = oi + bit' 1 . At an earlier epoch Z2 the crit- 
ical density is represented by barrier B2 = a 2 + &2t' 7 - For 
earlier look-back epochs or redshifts (increasing z), 02 > a%- 
The precise dependence of a on z is determined from the 
cosmological model. For example, in an fl m = 1 universe 
ffli,2 = 1.686 x (1 + Z1.2) [Note: the parameter a should not 
be confused with the scale factor which is often represented 
by the same notation] . 

The jagged curve shows a random walk trajectory that 
meets the barrier B\ for the first time at some scale t\ and 
the barrier B2 for the first time at some 'smaller' scale t2- 
This represents a collapsed object of scale t\ at redshift z\ 
that had a collapsed (or 'formed') progenitor of scale t2 at 
redshift 22. Let /(ti PI t2) be the probability of this occur- 
ring. To determine f(t\ n £2), we need to know the prob- 
ability that a random walk starting at the point O reaches 
the barrier B2 for the first time in time t2 — ti. In the new 
coordinates x,t, the barrier takes the form a + b(t + ti) 7 , 



where a — (12 — a\ — bit 1 and b = 62. Thus, 
f(ti n t 2 ) — g(a 1 ,b 1 ,0;t 1 ) g(a 2 — 01 — Mi, &2,*i; *a 
The conditional mass distribution f(t2\ti) is then given by 
f(t 2 I ti) = g(Aa - bxtj, b 2 , h; *2 - h), (8) 



(7) 



where Aa = 02 — ffli- Thus, a knowledge of g(a, b, to; t) solves 
both the conditional as well as unconditional mass distribu- 
tions. This was the motivation for introducing the parame- 
ter to into the modified barrier. In what follows, we will be 
considering the case when 61 = 62 = b. 

The first step in solving equation |5] is to make the 
boundary condition (equation[SJ time independent by choos- 
ing suitable coordinates. Let x — > B\(t) — x, with t remaining 
unchanged. Then, equation |S] reduces to 



dP{x,t) 



+ b~/(t + to) 



7-1 



dP(x,t) Dd 2 P{x,t) 



dt dx 
satisfying the boundary conditions 

P(0,t) = 0, 
P(oo,t) = 0, 
P(x,0) = Sd[x- (a + bt 1 )}, 



dx 2 



= 0, (9) 



(10) 

(11) 
(12) 



where So is the Dirac delta function. The first crossing dis- 
tribution is 



g(a,b,t ;t) 



D 
2 



d_ 

dt Jo 
dP 



P(x, t) dx. 



dx 



(13) 
(14) 



where the integral in equation 1131 has been evaluated using 
equation [5] 

We now bring equation |5| to a more convenient form by 
the following transformations. Let P(x,t) = 4>{x,t)e h( - x ' t \ 
Substituting into equation we determine the function 
h(x,t) by eliminating the terms proportional to d(f)/dx. 
Thus, if we choose h(x,t) to be 



h(x, t) = 



b^x{t + t ) 



7-1 



D 2£>(27 - 1) 

then equation |5| simplifies to 



d(j> 67(7 — l)z 
~dt + D{t + t ) 2 -~i' 



2 dx 2 ' 7 ^ 2 



(15) 



(16) 



When 7 = 2, the second term in equation llfcil becomes 
independent of time, and the equation becomes separable. 
We present the full solution for this case in section [1] Also, 
it turns out that when 7 = 1/2, the equation |U] can be made 
separable by a different set of transformations. We present 
the complete solution for this case in section|S| For arbitrary 
7, it is not possible to transform eouation ll6l into a separable 
form. However, it is possible to derive some results in the 
limit of large time. We discuss this in section [2] 



3 THE SQUARE ROOT BARRIER 

In this section we solve for the first crossing probability for 
the square root barrier (7 = 1/2). 



© 0000 RAS, MNRAS 000, 000-000 



4 A. Mahmood and R. Rajesh 



3.1 Solution 

Let y = xj y/t + to and 77 : 



. dP b- 



t. Then equation [5] simplifies to 
-ydP D d 2 P 



= 0. 



(17) 



dy 2 dy 2 

Equation 1171 is now separable. Let P(y,rj) = 0(77)^(2/) 
Then, 

(77 + 1 ) dQ 



d a 2 * (y 



efa) 977 
6) a* 



= -A (18) 

2* dy 2 ' 2* 9y +A ' (19) 

where A is an eigenvalue. The first of these equations is easily 
integrated to give 

1 



9(77) 



(20) 



(77 + io) A ' 

In equation 1191 we substitute (y — b)/v~D = f. Then, 
*(C)" + C*(C)' + 2A*(C) =0. (21) 
Further, letting 2A — 1 = v, and applying a transformation 
*(C) = *(C) ex P(-C 2 / 4 ), we obtain 

1 1 



* (0 + («+ 3 - 5 O*(0 



0. 



(22) 



The solutions to this equation are the parabolic cylinder 
functions U v (() and the eigen values v are determined 
by the boundary condition U v (—b/v~D) = (from equa- 
tion For large £, the function U V (Q goes to zero as 
C exp(-C 2 /4)[l + 0(l/()], consistent with the other bound- 
ary condition. 

It is interesting to note a similarity between the present 
problem and the problem of a quantum particle trapped in 
a potential V(Q given as 



no 



7S> 



ic 2 c > b — 



(23) 



In the limit 6 = 0, the problem reduces to the harmonic os- 
cillator potential with an V(£) = 00 for £ 0. The parabolic 
cylinder functions then reduce to Hermite polynomials with 
U 2n +i(z) = 2- n - 1 / 2 H 2n+1 {z)e-* 2 /\ 

For finite 6 exact eigen-values can be obtained nu- 
merically. Alternatively, the large n behaviour of v n may 
be obtained from the large v asymptotic behaviour of the 
parabolic cylinder function (see appendix 0. Then, 



2n- 



2y/2b 
ny/D 



yfn- 



2b 2 
'% 2 D 



+ (2fe-l) + 



(24) 

where k is an integer. 

The probability P(x,t) can now be written as a linear 
combination of all eigen modes as 



P(x,t) = £ 



(t + i )(»»+i)/a 



X — by/t + to 
y/D(t + to) 



X6XP [ 4D(t + t ) \- (25) 

Here A Vn are constants that are determined through the ini- 
tial condition. Using the orthogonality relation for parabolic 
cylinder functions 



U»(()U v (Od( = 0, a* j£ v. 



(26) 



10"' 

1 cr 2 
2 io~ 3 

1 cr 1 

10~ 2 
1 o~ J 





f 5 c = 2.0 


5 C = 3.0 


<5 C = 4.0 



10~ 2 10 1 10 10~ 2 10 1 10 

v = <5 C 2 /S = a'/t 

Figure 2. The unconditional mass function or the first crossing 
distribution for 7 = 1/2 (solid line) is compared with results from 
Monte Carlo simulations (histogram) for different values of S c . 
The dashed curve is the mass function obtained from the Press- 
Schechter constant barrier model. 



and using the initial condition 1121 one obtains 

2 /(4Di ) 



■u u 



A ^ Vr> j M (-6/vs) 

where we have defined I^—bjyS) as 



'Dio J 



/„ I -J=) = [°° uliQ d(. 

>DJ J-b/s/D 



(27) 



(28) 



Having obtained P(x,t), we can derive the distribution 
g(a,b,to;t) in a straightforward manner using equation 1141 



g(a,b,t Q ;t) = 



E 



to 



-b 2 /(iD) e a 2 /(4.Dt ) 

2(t + t ) X 
v/2 U' v {-b/y/D)U v {a/y/DTo) 



t + t 



I v (-b/VD) 



(29) 



where U' v (x) denotes a first derivative with respect to x. The 
conditional first crossing distribution can be obtained by 
using equation |H] To obtain the unconditional distribution 
we take the limit to 



g(a,b,Q;i) 



exp( 



-> in equation 
-b 2 /-iD) 



2t 



E 



to obtain 

v/2 uu-b/ 



D) 



I v (-b/VD) 



(30) 

It is straightforward to show that one recovers the constant 
barrier answer from this expression by setting 6 = 0. 

In figure [5] we compare the analytical result for un- 
conditional mass function (equation 1301 with results from 
Monte-Carlo simulations. In terms of the barrier form given 
in equation^] the parameter values are 7 = 0.5, j3 = 0.5 and 
q — 1. In the four panels (i, ii, iii and iv), a = S sc (z) = 1, 2, 3 
and 4 respectively The solid curve is obtained by truncat- 
ing the series in equation 1301 to the first 30 terms. There 
is a marked difference between the moving barrier and the 
constant barrier especially for large times (or small masses) . 

We also compare the expression for the square root bar- 
rier with Monte Carlo simulations for a 7 = 0.6 barrier. 
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<5 C = 1 .0 


5 C = 2.0 


<5 C = 3.0 


<5 C = 4.0 : 


0") I 


(iv) I 


10~ 2 10"' 1 10 10~ 2 10~ 1 1 10 
v = 6 C Z /S = o 2 /t 



' 0" 

-o- 3 



Figure 3. Same as in figure |3 except that the histogram is now 
for 7 = 0.6. 



This is shown in ngure[!J]where the solid curve is the expres- 
sion in equation 1301 and the histogram is data for 7 = 0.6. 
The other parameters are same for both simulation and the 
curve. It is seen that the data lies close to the curve. In fact, 
if the parameters b and a were allowed to be played with, 
a much better fit can be obtained. This shows the difficulty 
in distinguishing between various 7's by looking at the mass 
distribution. 



3.2 Large time behaviour of the unconditional 
first crossing probability 

In this subsection, we consider the large time behaviour of 
the first crossing probability. In the large time limit, the 
leading order contribution arises from the first term in the 
infinite series in equation 1301 This corresponds to the vari- 
able v = a 2 1 (Dt) — > limit. The unconditional first crossing 
distribution behaves as f(t) ~ t~ B with 6 — 1 + vo/2, where 
no is the smallest eigenvalue. It is not possible to analyti- 
cally calculate vq for arbitrary 6. When b — > 00, the eigen- 
value problem reduces to the harmonic oscillator problem 
with lowest eigenvalue being 0. Thus one expects 9 — > 1 as 
b — > 00. When b — > — 00, the barrier crosses the walker al- 
most immediately. Thus, one expects that 9 — ► 00 when 
— * —00. For small b (b 2 <C D), using the variational 
method, iKrapivskv fc Rednerl il996fl obtained an approxi- 
mate value of the exponent as 9 m 3/2 — b/v2irD. At the 
same time, one could naively assume that the large v be- 
haviour for parabolic cylinder functions (see eouation lA2ll is 
valid for all v and then estimate vo- Under this assumption, 
one obtains 



-7T 2 L> + 



7Y 2 D 



(31) 



Surprisingly, this expression agrees very well with the real 
answer and also has the right limits for large absolute values 
of b. Figure 0] plots the behaviour of the exponent 9 as a 
function of b/yD. This scaling is markedly different from 
the fixed barrier scaling in large time limit f(t) oc t~ 3 ' 2 . 

We mention here that, if we were to look at the large 
time behaviour of f(t) in equation with 7 = 1/2, then 



-6 -4 





b/VD 



Figure 4. The variation of the exponent 9 with b/\/T) is shown for 
the square root barrier. The dotted line is obtained by numerically 
solving U v ( —b/^/~D) = for v; the solid line shows the exponent 
obtained bv IKrapivskv fc Rednerl il996l) and the dashed line is 
from equation 1311 As expected, when b/y/D = 0, the value of the 
exponent is 3/2. 



9 = 3/2 — 7 = 1. Clearly, the lack of dependence on the 
parameter 6 makes it qualitatively different from the actual 
answer. 



3.3 Small time behaviour of the unconditional 
first crossing probability 

In the small time limit [v = a 2 /(Dt) 3> 1], we can deter- 
mine the asymptotic behaviour of the first crossing proba- 
bility g(a,b,0;t) by the saddle point method. In this limit, 
the leading contribution to g(a, b, 0; i) comes from a term 
with large v in the infinite series of equation 1301 On sub- 
stituting the large v asymptotic form of functions U' v and 
/„ (see Appendix A), and replacing the summation by an 
integral over n by noting from equation IA3I that dn = 
[1/2 + b/(2iv\/vD)]dv, equation 1301 reduces to 



/(*) = 



exp(-6 2 / 4I) ) 

T 



dv 



a 

Dt 



,/2 



v/2 



[1 + 0{1/v) 



(32) 



Here we have only retained terms up to 0(l/y/v). Rewriting 
the 'sine' term as an exponential, and evaluating the integral 
by the saddle point method, we obtain, after some algebra, 
the small time asymptotic form of the expression as 



/(*) 



■ exp 



(a + bVt) 2 



V2tTD^/ 2 

b(b 2 + 6D)Vt 
+ 24Da 



2Dt 



(33) 



Comparing this expression with equation it appears 
that the latter should be considered as an expansion for the 
small time behaviour of fit). Specialising to 7 = 1/2, f(t) 
has the same form as in eouation l33l except for the coefficient 
in front of the correction term. But whether this qualitative 
similarity continues to exist for 7 > 1/2 is not very clear. In 
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Figure 5. The analytical conditional mass function (solid curve) 
is compared with the results from Monte-Carlos simulations (his- 
tograms). The dashed curve is the fixed barrier conditional mass 
function. In panel (ii), the two curves are not distinguishable from 
each other. Here M\ is such that t\ = S(Mi). 
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Figure 6. The conditional mass function plotted in terms of 
variable v = Aa 2 /(t 2 — t\). The histograms show the results of 
Monte-Carlos simulations. The dotted curve shows the distribu- 
tion obtained by re-scaling of the unconditional mass function. 
The dashed curve is the constant barrier mass function. 



fact, in sections and we will show that equation [5] does 
not have the right form for arbitrary 7. 



3.4 Asymptotic analysis of the conditional mass 
function 

In this subsection, we study the asymptotic behaviour of the 
conditional mass distribution. This is given by equation |H] 
Using the result for g(a, b, to; t) (see eauation l29H . we obtain 



/(*a I *i) = T7T ex P 



bAa 



E 

{v} 



2t 2 

tiy /2 ui( 
t 2 



(Aa) 2 

4Dti 2D^/tl 

b/VD)U v [(Aa 



I v (-b/VD) 



(34) 



where Aa = a% — a-y. Unlike the fixed barrier case, the con- 
ditional mass function is no longer a universal function of 
just one scaling variable v — Aa/y/fo — £1. 

In figure El the conditional mass function as obtained 
from equation 1341 (solid curve) , is plotted with the results 
from Monte-Carlo simulations (histograms) for a range of 
present day halo masses (Mi) and look-back redshifts. Here 
M* is such that S(Af») = (5 2 c (0), which, for ACDM cosmol- 
ogy used here, has a value M* ~ 2 x lO 13 Af . The plot 
confirms the correctness of our expression. The result from 
the constant barrier model is shown for comparison (dashed 
curve). In figure|!)]we show the effect of rescaling the condi- 
tional mass function in terms of variable v. As expected, the 
shape of the curve is different at different look-back epochs 
and for different present day halo masses. 

We will consider the following limiting cases. Consider 
the limit (Aa/yDti) — * 00. This is the same limit as the 
one that was taken to obtain the unconditional mass distri- 
bution. Doing so, one obtains 



lim f{t 2 \t x ) = 



-b 2 /(4D) 



^ \Dt 2 



v/2 



UU-b/VD) 
Iv(-b/ 



D)' 
(35) 



which is the same as equation I3UI with a replaced by Aa. 
Therefore, for a given look back epoch, the conditional 
mass function tends to the unconditional mass function for 
large present epoch masses. Alternatively, for a fixed present 
epoch mass, the conditional distribution tends to the uncon- 
ditional distribution at large look back epochs. 

We now consider the opposite limit (Aa/y/Dti) — > 
0. Taking the Taylor series expansion of the function 
U(Aa/\/Dti — b/vD) we can write 



Aa - bJt 



U v {-b/VD) + -^=U' v {-b/VD), (36) 



where have only taken terms up to order Aa/VSti. Note 
that the first term on the right hand side is zero from the 
boundary condition. Thus equation 1341 can be written as 



f(t2 I h) 



Aa 



■E 



' J/2 UU-b/VD) 2 



I v {-b/VD) 



(37) 



In this limit there does not appear to be a direct way of dis- 
cerning the behaviour of /(ta However, figure||j]suggests 
that the form of the conditional mass function in this limit 
is closer to the constant barrier mass function. This point 
has also been discussed in ST02: for a given ti, at small 
look back epochs most random walkers will reach the second 
barrier in a relatively short time so that the barrier hasn't 
moved much and is effectively fixed. ST02 argue that the ef- 
fect should be more pronounced for massive haloes. Instead, 
panels (i) and (iii) of figure |H] suggest that for small mass 
haloes the conditional distribution are more like the constant 
barrier distribution, whereas for more massive haloes, the 
conditional distribution tends to the unconditional distribu- 
tion. Within the context of the present analysis, this may 
be understood by noting that for a given Aa, as Mi — > 00 
or equivalently t\ — > 0, equation 1341 approaches the uncon- 
ditional mass function. A simple argument can also be put 
forward by looking at the barrier form under consideration: 
B = Aa — b(tj — tj) [ see section|!2]. When t 2 — > ii (note that 
the minimum value of t 2 is ti) the situation is similar to 
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the fixed barrier case. On the other hand when t2 2> ti, the 
situation is similar to the single barrier problem. Therefore, 
for a given present halo mass and look back redshift, one 
expects a transition from a fixed barrier distribution to the 
unconditional moving barrier distribution. 



4 THE QUADRATIC BARRIER 

In this section, we solve for the first passage distribution for 
the quadratic barrier (7 = 2). We start with equation 1161 
which when specialised to 7 = 2 is 



2bx 



2 dx 2 



Let 



26 



(x,t) = 4>(x)T(t). Then, 

IdT 
T~dt ' 
D d 2 ip 
' 2^'dx 2 



-A, 



+A, 



(38) 

(39) 
(40) 



where A is an eigenvalue. 

The time dependent part has a solution T(t) = e" 
Let y = k(2bx/D — A), where 

k = [D/(2b 2 )] 1/3 . 

Then equation 1401 reduces to the Airy equation 

dy 2 



^2 -v$ 



: 0. 



(41) 



(42) 



The solution to the above equation are Airy functions Ai(y) 
such that when y — > 00, Ai(y) — > 0. The other boundary 
condition ij){x = 0) = implies Ai{— fcA) = and this fixes 
the eigen values A„. Hence, we can write solution in the form 



(x,t) = 22 C„e 



'Ai 



f2bx 



(43) 



where the coefficients C„ are to be determined from the 
initial condition. 

The function P(x, t) can now be written as 



P(x,t) = y^C n exp 



-X„t+ 



2bxt 2b 2 r 



D 3D 
The initial condition now gives us 
S D (x - a) = ^2 C n Ai 



A, 



A, 



(44) 



(45) 




Figure 7. The unconditional first crossing probability for a 
quadratic barrier. The fiducial values of parameters are b = 0.2 
and a = 1.675. The solid curve is from cquation l48l and histograms 
are from Monte-Carlo simulation of random walks. 



Now using equation 1141 we can write the expression for 
the first crossing distribution as 



/(*) 



2fe 2 fc 2 e- 262t3 /( 3 °> 
D 



E 



e . Xnt Ai{2ba/D-\ n ) 



Ai'(-kXn) 



Furthermore, the asymptotic form of the function Ai(— x) 
for x — > 00 is given as 



Ai(—x) 



,1/4 ' 



^ 8/2 , t 



(49) 



Thus using Ai{— k\ n ) = 0, we obtain the eigen values as 



3tt 
2 



n ~4 



2/3 



(50) 



In the large time limit t — > 00 it can be seen that 
the leading order term in the first crossing distribution is 
fit) oc exp[—2b 2 t 3 /(3D)]. This can be compared with the 
expression in equation [5] f(t) oc exp[— b 2 t 3 /{2D)}. Thus, 
when 7 = 2, the expression by ST02 is wrong. Figure |7| 
shows the first crossing distribution for 7 = 2 for fiducial 
parameter values: b — 0.2 and a = 1.675. The solid curve 
shows the result from eo uation 1481 and the histograms show 
the distribution obtained from Monte-Carlo simulations. 



Multiplying both sides of this equation by Ai [fc (2^2 — Am)] 
and integrating between the boundaries x = and x — 00, 
we obtain, using the orthogonality of Airy functions, 



At 



fc(^-A„ 



— Crr 



dxAi 



fc(^-A„ 



(46) 

The integral on the right hand side can be evaluated 
using the relation Ai"(x) = xAi(x) and is equal to: 
2^ Ai' 2 (— fcAm), where Ai' is the derivative of the Airy func- 
tion. Thus we obtain the coefficients as 



On 



2bk 



1 



D Ai' 2 (-k\ m ) 



Ai 



(2ba 



(47) 



5 ARBITRARY 7 

For 7 different from 0, 1/2, 1,2, it is not possible to obtain 
the full solution for f(t). However, it is possible to analyse 
it in the limit of large t. This analysis would be useful to 
test whether equation [5] has the right form or not. In this 
section we estimate 



lim 



Hf(t) 

t 2j-l > 



(51) 



such that fit) oc exp(crf 27-1 ). For comparison, equation [5] 
gives a = —b 2 /2. 

We start with equation 1161 with to — 0. cj>(x,t) obeys 
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the boundary condition <j>(0, t) 
bution f(t) is then given by 



0. The first passage distri- 



fit) " - ^ 
tW 2 dx 



exp 



(52) 



2D(27 - 1 

We will argue that the contribution from d(f> / 'dx\ x -o is at an 
order much smaller than the term in the exponential in the 
limit t — > oo. 

For large t, one could treat the potential in equation 
1161 as a slowly varying linear potential. We then make the 
adiabatic approximation. As t — > oo, we expect the system 
to be in ground state of the time dependent linear potential. 
Then, <j)(x, t) would be approximately equal to the first term 
in equation 1431 with b replaced by — l)t 7-2 /2. Thus, 
one would expect that the contribution from <j> to f(t) is 
utmost of the order exp[r 27_1 " 3 ], From equation 1521 we 
then conclude that 

i2 2 

— o 7 



" " 2,0(27 - 1) ■ (53) 

We now compare the above heuristic result for a with 
results from numerical simulations (with D — 1). It turns 
out that it is more convenient to numerically measure the 
survival probability S(t) rather than directly measure f(t). 
S(t) is the probability that a random walker has not crossed 
the barrier up to time t and is related to f(t) by f(t) = 
—dS(t)/dt; hence it has the same a as f(t). To estimate a 
we need to go to large times, but then face the problem that 
S(t) becomes exponentially small. It is difficult to overcome 
this problem using conventional Monte Carlo methods. In- 
stead, we use an algorithm kn own as "go with the winners" 
algorithm (iGrassbergerl 120020 . We briefly describe the algo- 
rithm. We start with N (N = 2 x 10 5 in our case) random 
walkers at the origin. When the number of walkers get re- 
duced by half (due to absorption at the boundary), copies 
are made of the remaining surviving ones, and the survival 
probability is halved. By repeating this procedure, we can go 
down to very low values of survival probability keeping the 
number of live walkers constant. As an example of the kind 
of data that is obtained, we show the survival probability 
for 7 = 0.75 in figure [S] 

From the numerical results for the survival probability, 
a is determined by best fit. In figure |5] we show the varia- 
tion of a with 7. The data agrees well with the formula in 
equation 1531 We would expect that equation 1531 is an exact 
result. 



6 HALO MASS FUNCTION FROM TV-BODY 
SIMULATIONS 

A two paramete r empirical formu l a for t he halo mass func- 
tion as given in ISheth fc Tormenl il999T> can be written as 



exp 



qv 
' 2 



(54) 



Here v = S BC (z) 2 /a 2 (m). The parameters are q = 0.707, 
p — 0.3 and A « 0.322, which is set by enforcing the condi- 
tion that integral of f(v) over all v equals unity. The Press 
Schechter formula has q = 1, p = and A = 1/2. For 
small v one obtains vf(v) oc u 5 ~ p . Equivalently for large 
t, f(t) oc * 3 / 2 - p . 
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Figure 8. The variation of the logarithm of the survival prob- 
ability S(t) with vt is shown. The dotted line has slope as in 
cquation l53l and is shown for comparison. Very small S(t) can be 
reached using specialised algorithms. 




Figure 9. The variation of a with 7 is shown. The values for 
7 = 1 and 7 = 2 are exact. The data is for b = 1.0. The solid line 
is the formula is equation 1531 



In this section we fit the TV-body mass function with 
the unconditional mass function for the square root bar- 
rier, taking b and q as free parameters (see equation^. For 
this, we use the small time asymptotic form of the uncon- 
ditional mass function as given in equation 1331 We use the 
halo catalogues of GIF simulations for ACDM cosmology, as 
available on http://www.mpa-garching.mpg.de/Virgo/ The 
cosmological parameters for the simulation are f2 m = 0.3, 
SIa = 0.7, h = 0.7 and as = 0.9. The box size of these 
simulations was L = 141 ft _1 Mpc and a total of 256 3 par- 
ticles were simulated. We only consider haloes with at least 
50 particles and a maximum mass of 3 x 10 14 M Q . This up- 
per limit is chosen as the statistics of haloes with masses 
> 3 x 1O 14 M0 is dominated by significant noise. For lower 
halo mass limit one should typically consider > 70 particles 
(Arif Babul private communication). The best-fit value of 
the parameters we obtain are b » 0.5 and q w 0.55. This 
value of q is somewhat smaller than that obtained by fit- 
ting the expression 1541 Interestingly, using the large time 
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v = (S c /a 



Figure 10. The halo mass function. Filled circles show the mass 
function obtained from GIF simulations. The error bars show the 
poissonian error. Left panel: The dashed curve shows the Press 
Schechter mass function and the solid curve shows the moving 
barrier mass function (see section||j]for details). The dotted curve 
shows the asymptotic form of the unconditional mass function 
( equation 1331 . Right panel: The s olid curve shows the b est-fit 
empirical mass function obtained in ISheth fc Tormenl il999l) . 

behaviour of f(t) as discussed in section 13.21 for b = 0.5 
one obtains f(t) oc t 3 ^ 2 ~ p , where p' « 0.2. In the following 
analysis we restrict to our best-fit values for q and b. 

In figure 1101 the left panel shows the GIF mass func- 
tion (filled circles) and the solid curve shows the mass func- 
tion as given in equation 1301 Also plotted for comparison is 
the Press Schechter mass function (dashed curve) and the 
asymptotic form of mass function as given by equation 1331 
(dotted curve). In the figure we have labelled the square root 
barrier mass function as 'ELL' as in ellipsoidal collapse. The 
fact that the full solution gives a worse fit than the asymp- 
totic formula (with the same parameters) shows the danger 
in fixing parameters through approximate formulae. 

The right panel shows the lSheth fc Tormenl Jl99Sl) mass 
function as given by equation 1541 It can be seen that the 
square root barrier mass function gives a reasonable descrip- 
tion of the halo mass function obtained from the simulations. 
It also shows the difficulty of extracting the correct 7 from 
the simulations. To do so, one needs to have good data for 
small mass, and the resolution at these mass scales is not 
good enough. 



7 IMPLICATIONS FOR MERGING HISTORY 
OF HALOES 

In this section, we present an application of our calculation. 
We use the results of the conditional mass distribution for 
the square root barrier to estimate the ma jor merger rates 
of dark matter haloes. In simple model sJWv ithe fc Loebl 
120021. 120031: iMahmood. Devriendt fc Silkll2004t) . halo merg- 
ers are followed by galaxy mergers. Major mergers of galax- 
ies are thought to trigger off quasars. Thus halo merger 
rates provide an estimate of quasar numbers at any given 
redshift. This is particularly useful, since quasars can be 
observed directly up to sufficiently high redshifts (z > 6). 
Previous studies have shown that the Press-Schechter (PS) 
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halo major merger rates are consistent with the evolution 
of qua sars at high redshifts (z > 2) but fail for low red- 
shifts JWvithe fc Loebll2002l) . The low redshift evolution of 
quasars is marked by an enhanced role of baryonic compo- 
nents and is therefore probably beyond the reach of simple 
predictions based on halo major merger rates. 

In the PS model the simplest estimates of halo 'cre- 
ation' rates are often taken as the positive term in the 
cosmological-time (or redshift) derivative of the uncondi- 
tional halo mass function (we explicitly write cosmological- 
time so as to avoid any confusion with the notation t, to, ti , 
used in the paper). This is partly a necessity as in the PS 
formalism all haloes are continuously accreting matter and 
are therefore newly formed. There is no direct way of differ- 
entiating between minor and major accretion episodes. To 
counter this deficiency of the model, a simple argument is 
used to define a major merger: the haloes which accrete mass 
comparable to their own mass, in a short interval of time are 
said to undergo a major merger. Using the conditional mass 
function, therefore, one can compute the probability that at 
a very small look back time Aa, a given halo had a 'formed' 
progenitor less than half of its present mass. Substituting Aa 
by da/dr (where r is the cosmological time), in this prob- 
ability yields the rate at which haloes of a given present 
day masses are forming through major mergers. We call 
these rates as major merger r ates of haloes. This approach 
has be en discussed, for e.g., in lMahmood. Devriendt fc Siikl 
12004) . 

Alternatively instead of looking back-wards one can 
compute the probability that a given halo merges with a 
comparable mass object ('for med' or otherwise) du ring a 
small step Aa forward in time JWvithe fc Loebll2003l) . The 
two approaches yield essentially similar results. 

7.1 'Creation' rates 

We first consider the implications of taking the positive term 
in the cosmic-time derivative of the unconditional mass func- 
tion. The left panel in figure ITTI shows the redshift evolution 
of creation rates as obtained from the Press Schechter mass 
function (dashed curve) and the mass function for the square 
root barrier (solid curve). These rates have been integrated 
over the minimum masses as depicted in the figure. It can 
be seen that the moving barrier mass function gives consid- 
erably higher creation rates for large redshifts. 

7.2 Major merger rates 

In order to compute major merger rates of haloes, we con- 
sider an object with a given ti (corresponding to mass Mi) 
at some epoch; for a small look back time Aa — + 0, we 
can apply the result of equation 1371 to obtain f(t2\ti). The 
fraction of objects with less than half mass progenitors can 
therefore be computed as 

p oc 

w(t 2 > t h \ t{] = / d* 2 /(t 2 |ti) 




where ti = S(Mi) and t h = 5(Mi/2). Replacing Aa by 
a = da/dr gives us w. This represents the fraction of ob- 
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Figure 11. Left Panel: The cumulative creation rates of haloes 
with minimum masses as depicted in the figure. Right panel: The 
cumulative major merger rates of haloes, above some given masses 
(see section [7. 21 for details). 



jects of mass Mi formed from major mergers, per unit time. 
Multiplying w with the existing number of objects of mass 
Mi as determined from the mass function, gives the 'ma- 
jor merger rate' of haloes of mass Mi. Integrating the rate 
over Mi from some minimum mass M m i n to 00 gives the 
cumulative major merger rate of objects above mass M m in- 
In the right panel of figure tTTl we plot cumulative major 
merger rates of haloes in the constant barrier (dashed curve) 
Press Schechter model and the moving barrier (solid curve) 
model. For the moving barrier model the unconditional mass 
function is as given by equation l33l Fieure lTTI shows that the 
moving barrier model yields a higher number of cumulative 
major merger rates towards high redshifts and particularly 
so for massive haloes (~ 10 12_13 Mq). This is precisely the 
halo mass range which is relevant for quasars. The inclusion 
of moving barrier threshold could, therefore, significantly 
affect the predictions of analytic and semi-analytic models 
dealing with the evolutionary history of the high redshift 
quasars. 

As an example we compute the number counts of 
quasars above a given redshift and flux level in observed soft 
(0.5 - 2 keV) and hard (2 - 1 keV) X-ray bands. For this, 
we us e the model discussed in lMahmood. Devriendt"fc" "sum 
<2004l) . The black hole mass (M B h) to halo mass (M h ) 
relation is as described in that paper. The major merger 
rates of haloes are computed as discussed above. Using the 
Mbh — Mh relation these rates are converted into black 
hole formation rates. The bolometric luminosities are esti- 
mated as the Eddington luminosity for a given black hole 
mass. To obtain the luminosity in observed X-ray bands, at 
a given redshift 2, we derive bolometric corrections in bands 
0.5(1 + z) - 2(1 + z) keV and 2(1 + z) - 10(1 + z) keV 
respectively. For this we use the spectral energy distribu- 
tion and the r e st-fram e bolometric corrections described in 
iMarconi etHI (12004) . The quasar life times are taken as 
the local dynamical times of galaxies (note that the rela- 
tive difference between the Press Schechter and the mov- 
ing barrier model will be independent of the lifetime). Thus 
luminosity functions in 0.5(1 + z) — 2(1 + z) keV and 
2(1 + z) — 10(1 + z) keV bands are obtained. 




Figure 12. Left panel: The number count of quasars in the ob- 
served hard X-ray band (2 — 10 keV). In each line style, the set 
of three curves from top to bottom correspond to z > 3, 5 and 
7 respectively. Solid curves show the moving barrier prediction 
and dashed curves show the prediction of Press Schechter model. 
Right panel: The number count of quasars in the observed soft 
X-ray band (0.5 - 2 keV). 

Quasar luminosities are related to the observed flux as 
Fx(Lx,z) — Lx /4:ir / Dl(z) 2 (here L x is the X-ray lumi- 
nosity). From the luminosity function <f)(L x ,z) of quasars, 
the number of X-ray quasars in the whole sky, above 
some flux level Fx an d above some redshift z is given as 
Jrlaiman fc Loebll 19981) 

f°° f°° ( d 2 V \ 

N(>F x ,>z)=4tt dz' dL x [——)<t>{L Xl z). 

Jz Jl(f x ,z) \dzd\lj 

(56) 

Here d 2 V/dz/dQ is the comoving volume element per unit 
redshift and per unit solid angle. The left panel of figure H2l 
shows the number count in a 6' circle, for hard X-ray band. 
The right panel depicts the results for soft X-ray band. The 
solid curve is the moving barrier model prediction and the 
dashed curve is the Press Schechter prediction. It can be 
seen that there is a pronounced discrepancy between the 
two models at high redshifts. Even though the full nature of 
this discrepancy can only be known through a model based 
on Af-body merger trees, our results highlight the difference 
qualitatively. It may be added here that the difference in 
the models primarily owes to the difference in the halo mass 
function in the two models. 



8 SUMMARY AND DISCUSSION 

In this paper we studied the first crossing probability f(t) 
of one dimensional random walks across a barrier moving as 
t 1 . Complete analytical solutions for the square root barrier 
(7 = 1 /2) and the quadratic barrier (7 = 2) were presented. 
For arbitrary power law barriers, large time estimates of 
the first crossing probability were derived. We showed that 
the formula for }{t) as presented in lSheth fc TormerJ (I2002T) 
is not valid for general 7. We also presented a methodol- 
ogy for approaching the two barrier problem for conditional 
mass function of haloes. This is relevant for deriving the 
progenitor distribution of haloes. Unlike the stationary and 
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the linear barrier, the conditional first crossing probability 
does not follow a simple re-scaling of the unconditional first 
crossin g probability for other b arrier forms. 

In ISheth. Mo fc Tormerl (j200lTl . it was argued that 
7 « 0.6 is the barrier arising from ellipsoidal collapse. It 
was also shown the mass function from JV-body simulations 
could be fitted well by a barrier of this kind. In this paper, we 
fitted the iV-body data using the square root barrier. This 
barrier has the advantage of being analytically tractable and 
the conditional distributions are also fully known. As shown 
in section^ the numerical data is fitted well with the square 
root barrier too. Since we also derived large time behaviour 
of first passage distributions across barriers with arbitrary 7, 
one could ask the following question. Treating 7 as a free pa- 
rameter, can it be determined using the data from iV-body 
simulations? For small times (large masses), f(t) is domi- 
nated by the stationary barrier. It is only for large times 
that 7 plays a role. Thus, to determine 7, the iV-body simu- 
lations should have good resolution at small masses. This is 
currently not available, and it is not possible to determine 
7 from the simulations with any degree of confidence. 

We also compared the predictions for merger rates of 
haloes in the constant barrier model with that of the moving 
barrier model. The moving barrier model predicts a signifi- 
cantly higher merger rates for massive haloes towards high 
redshifts. This arises as a consequence of the fact that the 
moving barrier mass function yields more massive haloes at 
high redshifts, as compared to the standard Press Schechter 
mass function. In this regard an interesting exercise would 
be to compare the halo merger rates obtained from sim- 
ulations with the constant and the moving barrier merger 
rates. In the present context it appears that the prediction 
of the moving barrier mode l is consistent with the e arly for- 
mation of massive galaxies (Glaz ebrook et al]l2004) and the 
presence of high redshift quasars with massive host haloes 
iFan e t al. 2003). We also presented a calculation of X-ray 
quasar number count above a given flux level and a given 
redshift. In terms of relative abundances, we found that the 
moving barrier model predicts systematically more quasars 
towards higher redshift. At z > 6 the number reaches almost 
dou ble the standard Prcss-Schechter prediction. 

Ivan den Boschl (1200211 has shown that there is a differ- 
ence between the average mass accretion history of haloes in 
simulations a nd that derived from the sem i-analytic merging 
history trees llSomerville fc K olatt 1999J). In particular the 
A^-body simulations suggest an earlier formation of haloes 
than is inferred from semi-analytic trees. Earlier formation 
epochs are reminiscent of higher merging activity towards 
high redshifts, which is the case in the moving barrier model. 
In this context it is worth pointing out that our expressions 
for conditional mass function could be used to generate "im- 
provised" merging history trees. The form of equation 1341 
suggests that the task of drawing progenitors using the given 
expression may not be as easy as in the case of a constant 
barrier model. However, the analysis in section l3~4l indicates 
that when Aa/y/ti — > one can use a considerably sim- 
plified expression for the conditional distribution ffalti) 
(equation EZJ. Therefore, for a given present day halo with 
mass Mi (corresponding to ti) one can choose an appropri- 
ately small A a and use equation 1371 to draw the progenitor 
masses. Alternatively we have pointed out that the two op- 
posite limits for the conditional mass function are the con- 



stant barrier distribution and the distribution obtained by 
simple re-scaling of the unconditional mass function. Hence 
the merging history trees in the moving barrier model would 
be constrained by these limiting distributions. In a forth- 
coming work this issue will be investigated in further detail. 
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large v asymptotic form of the function U v (y) and obtain 
the final integral as 



I v (-b/VD)=n^(-) V 



■0(l/v) 



(A8) 



APPENDIX A: PARABOLIC CYLINDER 
FUNCTIONS 

In this appendix, we discuss some standard results per- 
taining to parabolic cylinder functions and their integrals. 
Parabolic cylinder functions are the solutions of the differ- 
ential equation An integr al representat ion of these func- 
tions for v ^ — 1 is given as terdelvilll953|) 

17, ( X ) = -\/^ 2/4 dt e-' 2/ V cos (xt - ^) . (Al) 
For large v, the asymptotic form of these functions is 
U v {x) = \ft {^) V/2 [cos (xy/v- ^) 

+ 1^K-f) +o H- (A2) 

For the problem discussed in the paper, the boundary con- 
dition is U v (—b/\/~D) = 0. For large v, we can keep just the 
first term in eauation lA2l and estimate the eigenvalues to be 
given by 

= + $ -(* + .)!, (A3, 



where n is a large positive integer. Solving, we obtain 
v n = 2n-^:V^+^ + (2k-l) + o(^-]. (A4) 

Here, k is an arbitrary integer. 

To compute the asymptotic behavior of U' v (x) we note 
the recursion relation 2U' v (x) = xU v {x) — 2U v +i{x). Thus 
for x — —b/yD we have U' v (—b/v~D) = —U v +i(x). Using 
equation IA2l we can then write the large v asymptotic form 
for U' v {-b/sfD) as 

«U)-^(!r*(T+$>+«w)]. 



(A5) 

In order to compute the large v as ymptotic form of the 
integral I v (—b/V~D) we first note that iErdelvilll95'3T) 



Ut{y)dy= ^ 



V>(l/2 - v/2) - ip(-v/2) 



r(-v) 



(A6) 



where T(x) is a gamma function and ^>(x) is a logarith- 
mic derivative of the gamma function, defined as ip(x) = 
r'(a;)/r(x). Using a set of standard result s [for e.g. given 
in the appendix of iBender fc Orszad il978l) l we obtain the 
large v behavior of this integral as 

J U 2 v (y)dy = n^e{^) [l + 0(l/v)]. (A7) 

Now in the remaining integral J°^ b ^^j^U 2 (y)dy we use the 
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